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Abstract. - Hydrogen bonding is modeled in terms of virtual exchange of protons between 
water molecules. A simple lattice model is analyzed, using ideas and techniques from the theory 
of correlated electrons in metals. Reasonable parameters reproduce observed magnitudes and 
temperature dependence of the hydrophobic interaction between substitutional impurities and 
water within this lattice. 



Introduction. Hydrogen bonding [1-3] is responsible for the cohesion of water and 
is intimately involved in solvation. Protein secondary structures such as a helices and (3 
sheets are engendered by intramolecular H bonds that are stabilized in competition with 
water bonding. Such systems suggest that hydrogen bond correlations are relevant for their 
understanding. The purpose of this communication is to propose a simple toy lattice model [4] 
that allows consideration of proton correlations. We are motivated by the Hubbard model [5] 
for correlated electron systems, which has illuminated the behavior of transition metals, high 
temperature superconductors, heavy Fermion conductors, etc. Here we introduce for simplicity 
the most primitive version of our H-bond model, with molecular sites restricted to a simple 
cubic lattice. The phenomena of concern here involve short range effects, where the long range 
disorder of liquids is not paramount. 

The molecular structure of the water molecule is closely related to the unusually large 
electronic polarizability of the O ion. Pure coulomb considerations with a rigid oxygen 
atom would place the protons in a water molecule directly opposite one another, whereas the 
low energy configuration, due to the oxygen polarizability, locates them at a relative angle of 
104 degrees. We simulate this feature in a model which is simpler geometrically, with possible 
proton locations only at the six cartesian positions ±x, ±y, ±z of a simple cubic lattice. These 
positions, which we call basins, are located in the vicinity of the outer edge of the oxygen ion 
electronic cloud. The basins are meant to be local minima of the potential energy for the 
protons and would be determined by a balance between attraction to the nearest O and 
the crystalline electric fields associated with the nearby oxygens. We take the binding energy 
of a single proton in a basin associated with one of the O ion as the zero of potential energy. 
We treat the occupied basins as quantum states and the protons as spinless fermions, so that 
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any basin can be occupied by at most one proton. The potential energy cost for placing the 
second proton on an oxygen site already containing one proton will be V if the second site is 
perpendicular to the first (near neighbor), whereas the cost for putting the second opposite 
the first (next near neighbor) will be U > V. We consider a simple cubic lattice of such 
oxygen sites. Consider two neighboring water molecules (two protons on each) on that lattice 
in the x-direction (see Fig. 1). If there is a proton on the left molecule in the +x site, and 



Fig. 1 - Neighboring H2O molecules. The occupied proton basins are denoted by the small dark 
filled circles; the circle at the center on the right indicates the basin above the plane of the figure. 

the —x site on the right hand molecule is empty, then the proton can tunnel from the left 
to the right molecule, leaving a hydroxyl ion OHT on the left and creating a hydronium ion 
H 3 + on the right. The minimum energy cost for doing this is V, the value appropriate to the 
configuration illustrated in Fig. 1. We denote the tunneling matrix element by t. The notation 
has intentionally been chosen to draw attention to the close analogy of this with the well 
studied Hubbard model of correlated electron systems. The pH of water at room temperature 
suggests, by the law of mass action, a value of V ss 32fcsT, or about 0.8 eV. The virtual 
hop and return of a proton between neighboring molecules forms the hydrogen bond within 
the model, with a strength given in lowest order perturbation theory by —t 2 /V, suggesting 
t ps (5— 10)fcsT to yield the known strength of that bond at room temperature. The value of U 
will be greater than, but of the same order as V . Transitions, real or virtual, involving energetic 
cost U will be thermodynamically suppressed relative to those involving a cost V . Because 
it greatly simplifies the algebra without significant alteration to the physical conclusions, we 
will henceforth take U as effectively infinite, so that all corresponding configurations can 
be neglected. We note that we are also neglecting coulomb energies associated with ionic 
configurations. This may be of some consequence in the detailed physical behavior, but the 
number of real ions at temperatures of interest is extremely small. 

This picture also suggests a source of hydrophobic interactions [1-4, 6-9] . Consider replac- 
ing two of the water molecules by others without mobile protons, such as organic molecules. 
If these impurities are well separated, they each eliminate six hydrogen bonds, whereas if the 
impurities are near neighbors on the lattice, they eliminate only 11 such bonds. Energetics 
promotes clustering of the impurities, equivalent to a hydrophobic interaction. Of course, 
there are entropic effects, as well. Below we calculate the free energy, to second order in 
the hopping t, in order to estimate the size of this hydrophobic effect. The lattice model 
is, of course, most explicitly relevant to systems like clathrates or ice. The neglect of global 
reorganization of water in liquids overestimates the hydrophobic interaction there. 

We write the Hamiltonian H (including the chemical energy term —fiN) as the sum of a 
zcroth order term Hq, with no hopping (t = 0), plus the hopping perturbation H t . We will 
calculate the free energy from the log of the grand partition function, 




Z = 



Tre 



H 



(1) 



with a corresponding definition for the unperturbed grand partition function Zq. We have 
taken units of energy as fcgT here for simplicity of notation. The free Hamiltonian can be 
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written in terms of the parameters defined above as 

H o = jYl ( nia + n «*)( n n + n n) + ^ nia7lis ' ~ ^y^X n ic + n i&) ( 2 ) 

Here we have labeled the proton sites by 2 indices. The roman index i indicates the lattice site 
of the attached oxygen, and the greek index a labels the location on the oxygen: a = x, y, z 
and a — —a. The proton number operator on the site i, a is denoted by rii >a = c\ a Ci tOI . The 
hopping Hamiltonian is 



[4 +a , a c ha + H.c] . (3) 



A standard cumulant expansion gives as the second order correction to the free energy, 



Z - 1 



In— « / dX (H t (X')H t (X')) , (4) 

Z Jo Jo 

where we have used the standard notation for operators in the interaction picture, 

H t {X)=e XH °H t e- XH \ (5) 

and the angle brackets with subscript zero imply a thermal average determined by the unper- 
turbed Hamiltonian, 

(A) (Z )- 1 TrAe- H ° (6) 

Non-interacting free energy. - Since Ho is the sum over independent oxygen site Hamil- 
tonians, the corresponding partition function is the product of identical site terms for a 
lattice of N sites: Z — Zq . Moreover, we note that each of the terms in H depends only 
on N a = n a + n & , which simplifies the enumeration of terms. In principle, we should sum 
over proton number occupation on a site ranging from to 6, but the only terms which will 
play any important role are the obvious physical ones, namely those with numbers 1, 2, and 
3 (hydroxyl ion, water, and hydronium, respectively). Then we find 



6e^ + I2e 2f *- V + 8e 3tl - 3V (7) 



As always, the chemical potential n is determined by specifying the average number of protons, 
which here is 2 per site: 

2 = Am*, (8) 



which gives 



and the free energy 



Co«^/ 2 , (9) 



F /N = Ho-Tlnz . (10) 



Interaction energy; high temperature. - We must next calculate the second order pertur- 
bation correction to the free energy, given in Eq. ((4]). Again there are many equivalent terms. 
Each nonvanishing term in the trace (see Eq.©) involves the hop of a proton to a neighbor- 
ing oxygen site and its subsequent return to the initial site. Each of the cubic directions is 
equivalent, so we need to consider only a single pair of oxygens, with only a single pair of sites 
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between which the proton can hop. Every term in the trace depends on A and A' in the form 
exp[(A' — A)W], which gives for the integrals over those variables in (H|), 

1 r\ -W 1 1 

The energies W are positive or negative multiples of V, which is large compared to unity 
in magnitude at any temperature of interest to us. If W > 0, then the result is approxi- 
mately 1/W — 1/W 2 . If W < 0, then it is approximately e~ w /W 2 . If we neglect all terms 
exponentially small relative to the leading one, then we find the simple result, 

b | „ ™Vw«r> £ d ,£ m)H t^ ^ (1 - ±) . (i 2 , 

This represents the high temperature limit, in the sense that we require t 2 /V <C 1 (though 
we still assume V >■ 1), restricting the range of validity for physically realistic parameters to 
twice room temperature or more. We can include the leading corrections which depend on 
the chemical potential, which will be necessary to determine the free energy. The result for 
the partition function with all terms involving site proton occupancies of 1, 2, or 3 only, is: 

ln-^JV 1 [v W)+ 1 (13) 

Zo [6 + 3e-0>-v) + 4 e (M-2V)] 2 

With the fugacity given approximately by © we see that the limit (fT2"|) is given by the first 
factors in numerator and denominator here, with corrections of order e~ v / 2 multiplied by 1, 
V, or V 2 . At constant volume the (Helmholtz) free energy is given by F = Nfi — TlnZ, with 
the chemical potential /i determined by 

d 

N = 2N=—\xiZ. (14) 

In the correction terms of order t 2 we can take /i ~ ^jlq, as given by ([9]). Then we find for the 
fugacity Q = e M the approximate result 

C 2 - '^^e^. (15) 

We can immediately find the free energy (we re-insert all factors of the temperature T now 
for clarity) as 



F = 2Npi-TlnZ 



N 



V - —{1/V - T/V 2 ) - ±l e -W 2 - Tln(12) 
3 v3 



{fit < 1;/3V > 1). (16) 



The first two terms are the average energy up to second order corrections in the hopping, and 
the last term is the familiar form of the entropy in the high temperature limit for the set of 
12 states on each site where the two protons are near neighbors. Of course, we can trust this 
expansion only down to temperatures of order t 2 /V, well above the temperatures of physical 
interest for water. 

But in this high temperature regime we can readily estimate the hydrophobic interaction 
as described above, by looking at the correction to the free energy due to the hopping term t 
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in the Hamiltonian, which is what is eliminated locally when a water molecule is replaced by 
an organic or other molecule without mobile hydrogen atoms. We have 

SF = F - F = N(p - no) -T\n(Z/Z ) = NT\n((/( ) 2 - Tln(Z/Z ) (17) 

We have expressions for each of these terms, using (|9|), (fl5|) . and (fT3|) . If we keep only the 
leading exponential order of temperature dependence, we find 



3T 



V V 2 4 



(18) 



The hydrophobic interaction, then, is the difference in free energy between two impurities well 
separated and the same two as near neighbors, with one fewer virtual hopping opportunity 
excluded. We can read off this attractive interaction directly from (|18p : 



V 

V hydro ~ ^rp 



L + 11- y^_ p -0V/2 

V V 2 4 



(19) 



Low temperature limit. - We can also estimate the ground and low lying excited state 
energies to comparable accuracy to examine the low temperature limit and thereby to obtain 
an idea of the full temperature dependence of the interaction free energy which, as described 
above, suggests an origin of the hydrophobic interaction. 

In the absence of hopping, t = 0, it is clear that the ground state of the model consists 
of two protons in nearest neighbor sites on each oxygen. There is complete degeneracy with 
respect to the location of the independent pairs on each site. To second order in the hopping 
the most favorable situation corresponds to each hydrogen being able to tunnel virtually to 
a neighboring oxygen and then back, leaving in the virtual intermediate state a hydroxyl ion 
and neighboring hydronium ion, with the minimum possible extra intermediate state energy, 
namely V. In fact, we can construct such an arrangement (see Fig. 2). We take as a basis set 
for the doubly occupied sites the three near neighbor pairs xy, yz, and xz, which we denote 
0, +, and - for notational convenience. At the origin we place a molecule, at the position 
(1,0,0) a + molecule, and at (0,1,0) a - molecule. This configuration is repeated, with the 
basis vectors 110, 120, 011, giving three successive layers along the z-axis as shown in Fig. 2. 
Futher layers simply repeat this pattern. It is easy to check that this allows for every proton 
to participate in a tunneling process in which the intermediate hydronium has a proton in 
each of the three coordinate locations, so the intermediate energy is indeed V in every case. 
There is clearly some degeneracy, in addition to the obvious states arising from translational 
and rotational symmetry. That can be removed in higher order, but to second order in t this 
is the ground state energy of interaction, —2Nt 2 /V. 

We can further estimate the low temperature thermodynamics by examining the low lying 
excited states. There are the above-mentioned states which differ from the ground state only 
at higher order in the hopping t. But the number of these is not extensive (of order the number 
of sites N). There are, however, an extensive number of states in which a single site is occupied 
by a pair of near neighbor protons in such a way that the intermediate energy for hopping 
of one of them is U (which we have been taking infinite, or at least large enough to ignore) 
instead of V, and the interaction energy associated with the hopping of that proton is now 
—t 2 /U (zero, to be consistent with our earlier approximation for U) instead of — t 2 /V. There 
are 2A^ of these lowest excited states, obtained by taking a proton located in the direction 
w on an oxygen site in the ground state and moving it to the location —w, with w = x,y, 
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Fig. 2 - Ground state configuration in the z = plane. Open small circle lies below the plane of 
the figure. Successive planes above this one can be generated by shifting each row of the plane below 
down one row. 



or z. Thus there is an energy gap to the lowest (thermodynamically relevant) excited states. 
Since the low lying excited states are non-ionic (always with two protons on each site), we 
can calculate the free energy at low temperatures within the canonical ensemble: 



F = -ThiZ w N 



2 



V-\- Tlnfl -li, 1 



(/3t»l), (20) 



which decreases only exponentially slowly at low temperatures. Again, the hydrophobic re- 
pulsion corresponds to the blocking of one fewer virtual hopping sites when the impurities are 
near neighbors, so that interaction in this low temperature limit is 



V, 



hydro 



-2t 2 /y-Tln(l + 6e-' 3t2 / v ') (21) 



We plot the behavior of the hydrophobic interaction as a function of temperature in Fig. 3, 
using the limiting low and high temperature expressions above. It is convenient to scale both 
the interaction and temperature by the characteristic energy of the model t 2 /V: Vhydro = 
V hydro V/t 2 , and f = TV jt 2 . Then 



y= f2 + f]n(l + 6e-V*) (f«l); 
I (2/3T)(V/t) 2 (T»l). 



The result contains one factor of the parameter (V/t) 2 , which is of order 3 to 10. We have 
taken it as 6 for purposes of the plot, but the detailed result does depend on the specific value 
chosen. The results for the model do require an interpolation between the limiting curves, 
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Fig. 3 - The dimensionless hydrophobic attraction —V = Vh y dro(V/t 2 ) as a function of dimensionless 
temperature f = TV/t 2 , with {V/tf = 6 (see Eq. J22J). 

but we can conclude the hydrophobic interaction grows at low temperature to some ultimate 
maximum somewhat above room temperature, and then decreases. This suggests an increase 
in the hydrophobic interaction in the range of a few percent up to ten percent or so over 
a temperature change of, say 50 kelvin at temperatures where water is liquid, with values 
several times ksT, which is the correct sign and magnitude for what is observed. 

In summary, we have introduced a primitive Hubbard-like model to describe hydrogen 
bonding networks, that has the ability to rationalize the hydrophobic interaction in terms of 
disruption of this network. Natural extensions include hydration of ionic impurities, hydro- 
nium mobility, and competitive hydrogen bonding situations that occur in solubilization of 
polypeptides and polyethylene glycol. 
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